/* pde.h */
#ifndef EQPOISSON_H_
#define EQPOISSON_H_

#include <armadillo>

#define GAUSS_SEIDEL 		1
#define BLOCK_CYCLE_REDUCE 	2

typedef double (*ptr2func1arg)(double);				// real function of one variabe
typedef double (*ptr2func2arg)(double, double);		// real function of two variale
typedef ptr2func1arg f1v;							// short name for real function of one variabe
typedef ptr2func2arg f2v;							// short name for real function of two variabe

class EQPoisson
{
private:
	arma::vec BC_L;   		// container for discrete values of BC_left
	arma::vec BC_R;			// container for discrete values of BC_right
	arma::vec BC_T;			// container for discrete values of BC_top
	arma::vec BC_B;			// container for discrete values of BC_bottom
	arma::mat u;			// container for discrete values of u(x,y)
	double xlimits[2];		
	double ylimits[2];
	double eps_x;			// maximal acceptable step in x direction
	double eps_y;			// maximal acceptable step in y direction
	double hx;				// step in x direction
	double hy;				// step in x direction
	unsigned int Nx;		// number of points (length) in x direction
	unsigned int Ny;		// number of points (length) in y direction

	unsigned int _k;

	/* the Poisson equation is u_{xx} + u_{yy} = f(x,y), 0 <= x <= a, 0 <= y <= b; 
	 * with boundary conditions:
	 * u(x,0) = g_1(x) -------------> BC_left
	 * u(x,b) = g_2(x) -------------> BC_right
	 * u(0,y) = g_3(y) -------------> BC_top
	 * u(a,y) = g_4(y) -------------> BC_bottom
	 */
	
	f2v _f; 
	f1v _BC_left;
	f1v _BC_right;
	f1v _BC_top;
	f1v _BC_bottom;

	/* By the five-point difference method, it can obtain an algebraic equation 
	 *
	 * 							Au=b
	 *
	 * where A is a sysmtric banded matrix,
	 *
	 * 	   [ D   S             ]      	[ d   s             ]			[ 1   			]
	 *     [ S   D   S         ]		[ s   d   s         ]			[    1 			]
	 * A = [    ... ... ...    ],   D = [    ... ... ...    ],		S = [       ...		]/(h_y^2)	
	 * 	   [         S   D   S ]		[         s   d   s ]			[			1 	]
	 *     [             S   D ]		[             s   d ]			[			   1]
	 *
	 * and b = [b_1; b_2; ... b_{Nx-1}], d = -2(h_x^2 + h_y^2)/(h_x^2 h_y^2),  s = 1/(h_x^2)
	 */

	arma::mat D, S;
	arma::vec b;
	double d, s; 

	double determine_hx(int method);
	double determine_hy(int method);
	unsigned int compute_Nx(int method);
	unsigned int compute_Ny(int method);
	void init(int method = GAUSS_SEIDEL);
	void block_cycle_reduce_method();
	// double GS_iterate(int i, arma::vec x);
	void Gauss_Seidel_method(double precision);

public:

	/* parameter settors */
	void set_xlimits(double lower, double upper);
	void set_ylimits(double lower, double upper);
	bool set_f(f2v);
	bool set_BC_left(f1v);
	bool set_BC_right(f1v);
	bool set_BC_top(f1v);
	bool set_BC_bottom(f1v);
	void set_Nx(int);
	void set_Ny(int);
	void set_eps_x(double);
	void set_eps_y(double);
	void set_eps(double, double);
	bool solve(int method = GAUSS_SEIDEL);

	void save_data();
	void save_data(const char* filename) const;
	void draw_data();
	void draw_data(const char* filename) const;
};

#endif